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Abstract 


This  work  is  concerned  with  the  formulation,  implementation,  and  testing,  of  an  all  speed 
numerical  procedure  for  the  simulation  of  turbulent  mixing  and  evaporation  of  droplets.  The 
numeric  is  based  on  a  Multi-fluid  Eulerian  formulation,  with  a  pressure-based  fully 
conservative  Finite  Volume  method  equally  applicable  in  the  subsonic,  transonic,  and 
supersonic  regimes,  for  the  discrete  and  continuous  phases.  Three  two-equation  turbulence 
models  are  implemented  and  tested,  with  modifications  to  account  for  compressibility  at  high 
speeds.  The  models  are  the  k-e,  k-co,  and  more  recent  Shear  Stress  Transport  (SST)  models. 
The  code,  written  in  the  FORTRAN  computer  language,  is  supplied  with  a  number  of  ready 
to  run  cases.  Two  configurations  involving  stream-wise  and  cross-stream  spraying  are 
investigated  and  solutions  for  evaporation  and  mixing  in  the  subsonic  and  supersonic  regimes 
for  droplets  sprayed  in  laminar  and  turbulent  flow  streams,  are  generated  over  a  wide  range  of 
operating  conditions.  Results,  displayed  in  the  form  of  contour  plots  and  axial  profiles,  reveal 
the  degree  of  penetration  of  the  injected  droplet  within  the  supersonic  gas  phase,  and  the  rate 
of  evaporation  as  a  function  of  inlet  gas  temperature,  inlet  droplet  temperature,  and/or  the 
length  of  the  domain. 


Key-words:  evaporation,  mixing,  turbulence,  phase  change,  multi-phase  flow. 
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Nomenclature 

Aik) coefficients  in  the  discretized  equation  for  (f) k 1 . 

Bp  1  source  term  in  the  discretized  equation  for  (f>(k) . 

B(  k  1  body  force  per  unit  volume  of  fluid/phase  k. 

cfh,cfm  correction  coefficient  for  mass  and  heat  transport  in  droplet  model. 
C(p  1  coefficient  equals  to  1/  Rik)T(k) . 

CD  drag  coefficient. 

Ci,C2,C3  turbulence  model  constants. 

D (p)\(/>{k)]  the  Matrix  D  operator. 

Fd  drag  force. 

Fi,F2  blending  functions  used  in  the  SST  turbulence  model. 

Hp[ (f>{  k  ]  ]  the  H  operator. 

HP[uw]  the  vector  form  of  the  H  operator. 

I( ' )  inter-phase  momentum  transfer. 

j(pD  diffusion  flux  of  across  cell  face  ‘f . 

j'/l,c  convection  flux  of  <pik)  across  cell  face  ‘f . 

md  mass  rate  of  droplet  evaporation. 

Md  volumetric  mass  rate  of  droplet  evaporation. 

MWk  Molecular  weight  of  fluid/phase  k. 

P  pressure. 

Pr(jt),Pr  ^laminar  and  turbulent  Prandtl  number  of  fluid/phase  k. 
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q(k)  heat  generated  per  unit  volume  of  fluid/phase  k. 

Qik)  general  source  term  of  fluid/phase  k. 

Rlk)  gas  constant  for  fluid/phase  k. 

Red  Reynolds  number  based  on  the  droplet  diameter. 

Sf  surface  vector. 

Sc  Schmidt  Number, 

t  time. 

T  temperature  of  fluid/phase  k. 

U f  interface  flux  velocity  ■). 

u  velocity  vector. 

u,v  velocity  components  in  x-  and  y-direction,  respectively. 

X  Mole  fraction. 

Y  Mass  fraction. 

Greek  Symbols 

a  volume  fraction. 

J3(k)  thermal  expansion  coefficient  for  phase/fluid  k. 

8t  time  step. 

p  density. 

r  the  stress  tensor. 

A  conductivity  coefficient, 

r  diffusion  coefficient. 


<P 


general  scalar  quantity. 
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Ah  latent  heat. 

A p  [0(yt)]  the  A  operator. 

Ih^turb’^ff  laminar,  turbulent  and  effective  viscosity  of  fluid/phase  k. 
Q  cell  volume. 

Gi,02  turbulence  model  constants. 


Subscripts 


d 

refers  to 

E 

refers  to 

g 

refers  to 

k 

refers  to 

nb 

refers  to 

NB 

refers  to 

P 

refers  to 

s 

refers  to 

sat 

refers  to 

8 

refers  to 

CO 

refers  to 

vap,g 

refers  to 

the  droplet  discrete  liquid  phase, 
energy  equation, 
the  gas  phase. 

turbulent  kinetic  energy  equation. 

the  east,  west,  . . .  face  of  a  control  volume. 

the  East,  West,  . . .  neighbors  of  the  main  grid  point. 

the  P  grid  point. 

the  droplet  surface  condition. 

the  saturation  condition. 

turbulent  eddy  dissipation  equation. 

turbulent  eddy  frequency  equation. 

the  vapor  specie  in  the  gas  phase. 


Superscripts 
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C  refers  to  convection  contribution. 

D  refers  to  diffusion  contribution. 

(k)  refers  to  the  kth  fluid/phase. 

A  refers  to  variable  evaluated  at  reference  conditions. 

Old  refers  to  values  from  the  previous  time  step. 
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Introduction 


Background 

Recently  there  has  been  a  revived  interest  in  the  supersonic  injection  of  liquids, 
particularly  with  respect  to  fuel  injection  techniques  for  hypersonic  programs.  These 
designs  require  air-breathing  engines  capable  of  supersonic  combustion.  The  scramjet 
(supersonic  combustion  ramjet),  appears  at  present  to  be  the  only  practical  engine  for 
these  types  of  applications  since  it  is  able  to  produce  useful  thrust  at  hypersonic  flight 
Mach  numbers,  using  supersonic  flow  through  the  combustor.  Hydrocarbon-fueled 
scramjet  engines  are  of  particular  interest  because  of  their  dense,  non-cryogenic,  easily 
storable  fuels.  The  scramjet  concept  itself  is  fairly  old,  and  was  the  subject  of  studies 
throughout  the  1960s  and  again  in  the  1980s.  However,  its  coming  to  fruition  depends 
among  other  things  on  the  development  of  numerical  tools  for  the  simulation  of  its 
supersonic  combustion  process  and  related  phenomena.  More  specifically  effective 
disruption  and  mixing  of  hydrocarbon  fuels  in  a  supersonic  flow  is  a  crucial  ingredient 
for  the  success  of  any  scramjet  design.  Three  key  issues  govern  the  perfonnance  of 
the  liquid  injection  process  in  the  scramjet  engine,  namely:  the  penetration  of  the  fuel 
into  the  free-stream,  the  atomization  of  the  injected  fuel  drops,  and  the  level  of  fuel/air 
mixing  [1].  It  is  important  for  the  fuel  to  penetrate  effectively  into  the  free-stream  so 
that  the  combustion  process  produces  an  even  temperature  distribution.  If  the  fuel  does 
not  penetrate  sufficiently  into  the  free-stream,  the  combustion  will  mostly  occur  along 
the  surface  of  the  combustor  (or  injector  rake),  causing  inefficient  combustor 
operation  and  increased  cooling  problems.  Rapid  atomization  of  the  fuel  is  also 
required  for  efficient  combustion.  Increased  atomization  of  the  liquid  fuel  results  in 
increased  fuel/air  mixing  which  allows  a  higher  percentage  of  the  fuel  to  be  burnt  in 
the  short  time  before  the  entire  mixture  passes  out  of  the  combustor  (generally  on  the 
order  of  1  ms). 


Objectives 

Approaches  for  the  simulation  of  droplet  transport  and  evaporation  in  combustion 
systems  can  be  classified  under  two  categories,  namely  the  Lagrangian  and  Eulerian 
methods.  Within  both  methods  the  gaseous  phase  is  calculated  by  solving  the  Navier- 
Stokes  equations  with  a  standard  discretization  method  such  as  the  Finite  Volume 
Method.  In  the  Lagrangian  approach  [2,  3],  the  spray  is  represented  by  discrete 
droplets  which  are  advected  explicitly  through  the  computational  domain  while 
accounting  for  evaporation  and  other  phenomena.  Due  to  the  large  number  of  droplets 
in  a  spray,  each  discrete  computational  droplet  is  made  to  represent  a  number  of 
physical  droplets  averaging  their  characteristics.  The  equations  of  motion  of  each 
droplet  are  a  set  of  ordinary  differential  equations  (ODE)  which  are  solved  using  an 
ODE  solver,  a  numerical  procedure  different  from  that  of  the  continuous  phase.  To 
account  for  the  interaction  between  the  gaseous  phase  and  the  spray,  several  iterations 
of  alternating  solutions  of  the  gaseous  phase  and  the  spray  have  to  be  conducted. 
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Therefore,  the  computational  effort  for  strongly  interacting  two-phase  flows  with  the 
Lagrangian  method  is  rather  large.  Furthermore,  for  turbulent  flow  simulation  the 
above  model  has  to  be  augmented  with  a  stochastic  or  Monte-Carlo  approach.  In  the 
Eulerian  approach  [3,4],  the  evaporating  spray  is  treated  as  an  interacting  and 
interpenetrating  continuum,  in  analogy  to  the  continuum  approach  of  single  phase 
flows,  each  phase  is  described  by  a  set  of  transport  equations  for  mass,  momentum  and 
energy  extended  by  interfacial  exchange  terms.  This  description  allows  the  gaseous 
phase  and  the  spray  to  be  discretized  by  the  same  method,  and  therefore  to  be  solved 
by  the  same  numerical  procedure.  Because  of  the  presence  of  multiple  phases  a 
multiphase  algorithm  is  used  rather  than  a  single-phase  one. 

In  this  project  we  were  interested  in  developing  the  numeric  foundation  for  the 
simulation  of  supersonic  combustion,  more  specifically  to  develop  a  code  for  the 
simulation  of  supersonic  droplet  evaporation  and  mixing.  This  was  achieved  through  a 
multi-fluid  all  speed  pressure-based  finite  volume  flow  solver  in  which  a  droplet 
evaporation  model  was  implemented.  Three  two  equations  turbulence  models  have 
been  considered  to  account  for  the  droplet  and  gas  turbulence,  all  modified  for 
supersonic  flows.  These  are  the  k-e  [5],  the  k-co  [6]  and  the  SST  [7]  models  for  the 
gas  phase.  Droplet  turbulence  is  estimated  using  an  algebraic  model  based  on  the 
Boussinesq  approach.  The  model  as  it  stands  does  not  account  for  droplet  breakup  or 
coagulation,  but  different  droplet  sizes  are  accounted  for.  The  use  of  an  Eulerian 
approach  has  many  advantages:  same  validated  numeric  used  for  all  phases,  ease  of 
implementation  of  acceleration  techniques  (such  as  multi-grid),  improvements  to  code 
are  carried  over  to  all  phases. 

The  model  developed  and  implemented  in  this  project  accounts  for  the  heating  and 
evaporation  of  multiple  size  droplets.  The  interaction  of  the  gas  and  droplet  phases  is 
accounted  for  by  using  a  two-phase  multi-component  Eulerian  approach.  The 
implementation  has  been  done  for  multi-fluid  flow  at  all  speeds,  and  forms  a  solid 
foundation  for  the  future  study  of  the  more  general  problem  of  fuel  mixing  and 
combustion.  An  in-house  all-speed  multiphase  finite  volume  CFD  code  is  used  to 
compute  both  the  gas  flow  field  and  the  disperse  droplets  phase. 


Plan  Of  Report 


In  what  follows  detailed  derivations  of  the  governing  equations  for  droplet  transport 
and  evaporation  are  presented.  Then,  results  obtained  are  discussed.  Computations 
have  been  performed  for  subsonic,  supersonic,  laminar,  and  turbulent  flows. 
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Governing  Equations 

Interpenetrating  Media 


When  describing  multiphase  flow  phenomena  it  is  implied  that  more  than  one  phase 
exist  within  a  small  volume  at  any  particular  time.  This  view  rests  on  the  idea  of  time 
and  space  averaging  and  implies  the  following: 

1 .  Any  small  volume  can  be  regarded  as  containing  a  volume  fraction  a(k  1  of  the  kth 
phase,  so  that,  if  there  are  n  phases  then 

2>w= i  (i) 

k=  1 

2.  The  mass  flow  rate  of  any  phase  across  any  elementary  surface  within  the  domain 
at  any  time  can  be  expressed  in  terms  of  the  quantity  aik)  pik'fv{k)  where  pVh:)  is  the 

density  of  the  kth  phase  and  v(;)  is  the  local  value  of  the  velocity  vector  of  that 
phase. 

3.  When  the  content  of  the  finite  volume  and  the  flow  rates  across  finite  areas  are  to 
be  computed  over  finite  time  intervals,  a  suitable  averaging  over  space  and  time 
will  be  carried  out. 

Observance  of  these  rules  amounts  to  treating  each  phase  as  a  continuum  in  the  space 
under  consideration,  which  requires  choosing  a  proper  scale  for  the  control  volume 
used,  or  a  Representative  Elementary  Volume  (REV). 

For  a  multi-phase  system  the  equations  are  derived  over  a  REV  within  which  the 
different  phases  are  present.  The  scale  of  the  modeled  phases  has  to  be  smaller  than 
the  REV.  An  extensive  review  of  this  averaging  approach  can  be  found,  for  example, 
in  Hassanizadeh  [8,9]. 

Except  for  the  near  region  of  the  injector  where  the  spray  is  dense,  the  volume  fraction 
of  the  spray  is  low.  For  this  dilute  two-phase  flow  regime,  interaction  between 
droplets  can  be  neglected.  Starting  from  the  Navier-Stokes  equations,  instantaneous 
transport  equations  for  the  gas  and  droplet  phases  can  be  derived  either  by  spatial, 
temporal,  or  ensemble  averaging.  However,  these  transport  equations  can  only  be 
used  for  the  description  of  sprays  in  laminar  gas  flows.  Since  combustors  generally 
operate  in  the  turbulent  flow  regime,  the  system  of  equations  is  extended  by 
introducing  turbulent  fluctuations  of  the  transport  quantities  followed  by  Reynolds 
averaging  of  the  equations.  For  the  gaseous  phase,  three  two-equations  models  were 
employed,  namely  the  k-e  of  Spalding  and  Launder  [5,10],  the  k-co  of  Wilcox  [6]  and 
finally  the  SST  model  of  Menter  [7],  all  modified  to  account  for  compressibility 
effects  at  high  speed  flows.  An  algebraic  model  based  on  a  Boussinesq  [11]  approach 
approximates  the  turbulence  terms  in  the  droplet  phase  transport  equations.  The 
interacting  flow  fields  are  described  by  the  transport  equations  described  in  the 
subsequent  sections. 
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Geometric  Conservation  Equation 

The  volume  fractions  ag  and  0Cd,  for  the  gas  and  droplet  phases  respectively,  are  a 
result  of  the  averaging  process.  For  the  two-phase  flow  considered  they  are 
characterized  by  the  condition: 

a A  +  a%  =  1;  aA  «  1  (2) 


Gas  Equations 


Since  we  are  simulating  evaporation  in  the  turbulent  flow  regime,  the  system  of 
transport  equations  is  extended  by  introducing  turbulent  fluctuations  of  the  transport 
quantities  followed  by  Reynolds  averaging  of  the  equations.  The  continuity 
momentum  and  energy  equations  for  the  gas  phase,  composed  of  two  species  namely 
air  and  vapor,  are  generally  written  as  described  below. 


Mass  Conservation 


The  mass  conservation  equation  is  given  by: 

(«*  Ar )+ V-k  V>J  = =  Mvap,g  (3) 

where  ag,  vg,  and  pg  are,  respectively,  the  volume  fraction,  velocity,  and  density  of  the 
gas  phase.  The  density  is  computed  from  the  air  and  vapor  densities  as: 

J_ 

Pg  ~k= i 

where  Y  (/°  and  p  U  i  are  the  mass  fraction  and  density,  respectively,  of  the  kth  species, 
and  ~k=(air,vap)  indicates  a  summation  of  k  over  the  air  and  vapor  species. 

The  density  can  also  be  obtained  from  the  ideal  gas  relation  as: 


Y 


(*) 


(air, vap )  P 


(*) 


(4) 


Ps  = 


RgTg 


P 

[rj  mw g]rg 


Rn 


MW, 


f  Y(k)  ^ 


(5) 


Rr 


g  J 


y 

i  \  MW 

-k=[air  ,vap ) 


(*) 


In  equation  (3),  the  source  term  is  due  to  the  evaporation  of  the  fluid  droplets  and  is 
thus  equal  to  the  volumetric  mass  rate  of  droplet  evaporation,  thus  we  have: 


Mg=Mvap,g=-Md 


(6) 


Momentum  Conservation 


The  momentum  conservation  equation  of  the  gas  phase  is  written  as 

0  == 

^(wJ+V'kP,vgvg)=-^  +  V^+  FgB  +FgD  +  Mvapg\A 


(7) 
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where  FB  and  F°  are  body  and  drag  forces  respectively,  for  the  case  of  the  gas  phase 
the  body  force  can  be  neglected,  while  the  drag  force  due  to  liquid  droplets  is  written 
as: 


K IK  -  klk  ~vo) 


(8) 


4  s  D 

where  D  is  the  droplet  diameter  and  the  aerodynamic  drag  coefficient  is  given  by  [1 1]: 


_  .  . .  24  5.48 

CD  —  0.36  H - h 


Red  Ref73 


with  Red  being  the  droplet  Reynolds  Number  defined  as: 

pJk 


Re 


The  stress  tensor  is  written  as 

T  =  ugHeff  (Vv  +  Vvr  -f(V-v)l) 
where 

Aeff  =  A  +  M  tilth 


(9) 

(10) 

(ID 

(12) 


Turbulence  Models  in  the  Primary  Gas  phase 

The  modeling  of  turbulence  in  multiphase  flows  is  an  ongoing  research  area.  As  in 
single-phase  flows  there  is  no  particular  model  that  can  accurately  predict  turbulence 
over  a  wide  range  of  flows,  therefore  various  formulations  are  needed  to  cover  the 
numerous  types  of  multiphase  flows.  Furthermore  research  in  multiphase  turbulence 
is  still  in  its  infancy  as  compared  to  the  work  accomplished  in  single  flow  turbulence, 
and  turbulence  models  are,  understandably,  still  based  on  extensions  of  single  flow 
turbulence  models.  Of  the  different  groups  of  turbulence  models,  the  eddy  viscosity 
models  are  those  used  in  the  majority  of  turbulent  multiphase  simulations.  In  this 
respect,  two  different  strategies  are  employed:  (i)  a  dispersed  fonnulation  in  which  the 
turbulence  quantities  in  the  dispersed  phase  are  based  on  the  turbulence  model  of  the 
continuous  phase  [13-18];  and  (ii)  a  phasic  fonnulation  in  which  a  set  of  turbulence 
model  equations  are  solved  for  each  phase  [19]. 

As  the  focus  of  the  present  work  is  on  a  special  type  of  disperse  flows,  namely  the 
flow  of  sprayed  droplets  and  as  these  generally  operate  in  the  turbulent  flow  regime, 
the  system  of  transport  equations  need  to  be  extended  by  introducing  turbulent 
fluctuations  followed  by  some  averaging  procedure  (Favre  or  Reynolds  averaging). 
For  the  current  application,  three  two-equation  turbulence  models  (k-e,  k-co,  and  SST) 
are  employed  for  the  continuous  phase  to  represent  the  transport  tenns  resulting  from 
correlations  of  fluctuating  quantities,  while  for  the  dispersed  phase  the  turbulence 
tenns  are  approximated  by  an  algebraic  model  based  on  a  Boussinesq  approach  as 
described  later. 

Two-equation  turbulence  models  are  the  best-known  and  most  widely  used  models  in 
industrial  applications.  Amongst  these  the  k-e  model  of  Launder  [5,10]  is  the  most 
widely  used.  Although  a  number  of  additional  two-equation  models  have  been 
developed  over  the  years,  only  one  of  these  models  has  gained  significant  attention 
and  challenged  the  dominance  of  the  k-e  model,  namely  the  k-co  model  of  Wilcox  [6]. 
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For  both  models,  a  number  of  modifications  and  improvements  are  available,  which 
greatly  extend  their  applicability. 


k-£  model 


The  well-known  k-e  [5]  model  is  based  on  the  Boussinesq  approximation  with  the 
turbulent  viscosity  fonnulated  as: 


M, 


turb(k—£),g 


pCki 

rg'-'u  - 


(13) 


where  the  turbulent  kinetic  energy  kg  and  the  turbulent  energy  dissipation  rate  £g  are 
computed  using: 


d_ 

dt 

a_ 

dt 


(<*gPgkg)+v(agPgvgkg )  =  V.(agPeffAgVkg)+  ag 
(agPgeg  )+V.(agpgvs£g )  =  V  .{agPeff£  gV  eg )  +  ag 


(pk-pg£g)+sk„ 


(  E  E1  ^ 

CeX^Pk-C£lPgf 
v  j 


(14) 


+  5 


£,d 


where 


(15) 


Peff, 


k,g 


Ml 


am,g 


+ 


Mturb(k-e),g 


(16) 


Meff 


£,g 


Mlam.g 


Mturbik- 


■e),g 


with  the  following  model  constants: 


(17) 


T 


Cl 

c2 

c. 

Gk 

Gw 

Pr, 

Set 

1.44 

1.92 

0.09 

0.90 

1.30 

0.85 

0.7 

he  production  of  turbulent  energy  term  is  given  as: 


vVv 


(18) 


k-co  models 

The  k-co  model  of  Wilcox  [6]  is  similar  in  structure  to  the  k-£  model  and  is  also  based 
on  the  Boussinesq  approximation.  Two  transport  equations  are  solved  to  determine 
the  two  (large)  scales  of  turbulence.  The  advantage  of  replacing  the  £-equation  by  the 
co-equation  is  that  the  second  is  easier  to  integrate  (more  robust),  that  it  can  be 
integrated  through  the  sub-layer  without  the  need  for  additional  damping  functions, 
and  that  it  performs  better  for  flows  with  weak  adverse  pressure  gradient.  The 
conservation  equations  are  written  as: 


0 

jt  (agPgkg)+V-(agPgVskg )  =  v{agPeffAgVkg )+  ag  (pk  -  p*Pgkg(og)+  Sk4 


(19) 
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-  (agpgcog)+V.{ccgpgvgcog )  =  v{agpeff^gV  cog ) 

with  the  following  model  constants: 


f  CO  ^ 

Ca]-^Pk-Cfi]pgcog 

V  KS  J 


+  S. 


(oA 

(20) 


Cal 

Cpi 

Cn 

Oki 

Owl 

Pr, 

Sct 

5/9 

3/40 

1 

0.09 

2 

2 

0.85 

0.7 

where 


P turb(k - 


(0),g 


pg 


CO 


Peff,k,g  ~  Plain, g 


P 


turb{k-co),g 


'k\ 


Peff, 


CO,g 


'  Plam,g 


P 


turb(k-co),g 


co  1 


(21) 


The  major  drawback  of  the  Wilcox  model  is  that  it  is  very  sensitive  to  the  values 
specified  for  co  in  the  free  stream  [20],  which  leads  to  strong  dependence  of  the 
solution  on  the  arbitrary  specification  of  the  free  stream  co.  This  dependence  is  not 
present  in  the  k-e  model. 


The  SST  Model 

In  order  to  overcome  the  problem  of  free  stream  dependency  in  the  Wilcox  model, 
Menter  [7]  combined  the  two  models  so  as  to  take  advantage  of  their  respective 
strength  i.e.  the  robustness  of  the  k-co  model  near  wall  surfaces  due  to  its  simple  low 
Reynolds  number  fonnulation  and  its  ability  to  compute  flows  with  weak  adverse 
pressure  gradients  accurately,  and  the  better  performance  of  the  k-e  model  near  the 
boundary  layer  edge  and  away  from  walls,  due  to  its  insensitivity  to  the  free  stream 
values.  The  basis  of  this  technique  is  the  transformation  of  the  k-e  model  to  a  k-co 
formulation.  This  is  an  exact  transformation,  except  for  small  contributions  from  the 
diffusion  term  due  to  the  difference  in  the  diffusion  coefficients  of  the  k  and  the  e 
equations.  The  k-co  formulation  of  the  k-e  model  is  given  by: 


d_ 
dt 

Jt  (agPg<og)+v{ocgpgv%CQg )  =  V.for,// 


(ccgpgkg)+v{agpgvgkg )  =  V.{agpeffAgVkg )+  (pk  - (1  pgkgcog )+  S , 


+  cr 


g^eff,co,g  '  *"g 

Ca2CyLPk-Cp2pgco2g 
Ks  J 


+  2(7,,,-,  —  Vk  :  V  co„  +  S„ 


(22) 


(23) 


0)2  '  ,vg 

CO 


Jg  1  ^ co,d 


with  the  additional  constants  assigned  the  following  values: 


Ca2 

Cp2 

Ok2 

Ow2 

0.44 

0.0828 

1 

0.856 

The  BSL  model  is  derived  by  multiplying  the  k-co  equations  (Eqs.  (19), (20))  with  a 
blending  function  Fi  and  the  k-co  formulated  k-e  model  equations  (Eqs.  (22), (23))  by 
(1  -F i),  yielding  the  following  equations  for  k  and  co  [7]: 
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jt  k/?A )  +  V-kP,v A )  =  )+  ag (/>  - P*Pgkgcog)+  SkJ  (24) 

d 

-  (ccgPg(Qg )  +  V.{agPgvgcog )  =  V.(agPeff^gW  cog ) 


+  ar 


f  co 
a 

v  K 


s  p, 


Ppm 


2 

g~g 


+  29„  — Vt  -V® '(l-Fj+S, 


(25) 


co„ 


(0,d 


This  equation  is  fonnally  very  similar  to  that  of  the  standard  k-co  model,  however  all 
its  coefficients  depend  on  the  blending  function  Fi,  in  the  form: 


0  =  F101  +  (1  —  F>2  (26) 

The  blending  function  Fi  depends  on  the  solution  variables  and  on  the  distance  from 
the  nearest  wall  and  is  given  as: 

P]  =  tanhfy,4)  (27) 

where 


yl  =  Min 


Max 


#7  500v 

P  o)gy  y  ay 


CDkay2 


and 


(28) 


CDu„  =  Max\ 


d/cg  day 


g^col 


COg  dXj 


dx, 


,10 


-10 


(29) 


The  BSL  model  has  a  similar  performance  as  the  k-co  model  for  boundary  layer  flows 
and  is  nearly  identical  to  the  k-e  model  for  free  shear  flows.  Its  robustness  is  close  to 
that  of  the  k-co  model. 


A  further  modification  to  the  BSL  model  yields  the  SST  [7]  model  which  when 
compared  to  other  eddy-viscosity  models  has  an  improved  adverse  pressure  gradient 
performance.  The  modification  is  in  the  definition  of  the  turbulent  viscosity,  which 
takes  the  fonn: 


Pturb,  i 


(SST),g 


palkg 

Max[alO)g,Q.F2 ) 


_  .  Pturb,(SSTXg  Pturb,(SST),g 

Peff,k,g  ~  Plain, g  ~  Peff,m,g  ~  Plam,g  ~ 


where  Q  is  the  magnitude  of  the  vorticity,  and  F2  is  given  by  [7]: 


F2  =  tanh(y22 )  with  y2  =  Max 


2Jk^  500v 


n *  ’  2 

P  o)gy  y~co 


(30) 


(31) 


Conservation  of  Energy 

The  energy  conservation  equation  is  written  using  static  enthalpy  as  the  main  variable: 
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d_ 
dt 
with 


) + v-fc  p  syA ) : =  -v  •  %  +  v- 


rM!urtA_  Wh  ^ 

Pr,  g 


+  Sh,g  +  Mvap,gKap,s 


(32) 


(33) 

In  Eq.(33)  X  is  the  thennal  conductivity,  and  the  gas  enthalpy  is  defined  as: 

hg  =  HYk A  =  Yair^hairg  +  Yvapghvapg  (34) 

~k=(air,vap ) 

where 

hair,g  =  Cp,air,g  (Tg  ~  T  )+  Kir,g  (35) 

Kap,g  =  Cp,vap,g  ^g  ~Td)+  CP,dTd  +  Ah(Td)+  Kap,g 

The  Sh,g  tenn  on  the  RHS  of  Eq.  (32)  is  the  source  term  due  to  gas-droplet  interaction 
written  as: 

S^=nDY(Td-Tg)  (36) 

while  in  the  last  term: 

Kap,g,s  =  AKap  iTd  )  +  Cpjd  (37) 

with  the  subscript  “s”  being  used  to  indicate  that  the  enthalpy  is  evaluated  at  the 
droplet  surface  temperature,  which  in  the  case  of  the  Unifonn  Temperature  Model 
(adopted  in  this  work)  is  equal  to  the  droplet  temperature  Td.  Moreover,  T°  is  the 
reference  temperature,  h°  is  the  enthalpy  of  formation  at  the  reference  temperature,  Td 
is  the  droplet  temperature  (liquid  state).  Moreover,  it  is  to  be  noted  that  the  enthalpy  of 
the  vapor  contains  in  addition  to  the  enthalpy  of  formation  a  term  known  as  the  latent 
heat  of  vaporization,  A hvap(Td),  which  accounts  for  the  enthalpy  absorbed  during  the 
phase  change  from  liquid  to  vapor. 


Conservation  of  Species 


In  Eq.  (32),  Y  is  the  mass  fraction  of  vapor  in  the  gas  phase  and  is  obtained  by 
solving  the  following  mass  fraction  conservation  equation: 


ft^gPgYvapJ+V.(agPgygYvapg )  =  V.(agr%^VTra^g )+  agPgYvapg 
where 


(38) 


Ty 

*vaP’g  -<eff 


=  Y\ 


+  - 


M, 


turb,g 


(39) 


and  <7 Ya  is  the  turbulent  Schmidt  number  of  value  0.7  as  given  by  Spalding  [21]. 

Moreover,  the  source  tenn  on  the  RHS  of  equation  (38)  arises  in  the  conversion  of  the 
species  production  rate  to  a  concentration  rate,  as  shown 


Y, 


Yvap,g  (t  +  At)  -  Yvap  g  (t) 


vap,g ' 


vap,g 


At 


(40) 
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where 

Yvap,gW  = 


in,. 


m  +  m 

vap,g  air,g 


m  +  Am 

Yvap,g(t  +  At)  =  -  mp’8 


m  +  Am  +  m  . 

vap,g  vap,g  air,g 


and 
1  Y, 


Y. 


_  vap  air 


Pg  P  vap  P  air 

Substituting  equation  (41)  into  equation  (40),  we  get: 


m  +  Am 

vap,g  vap,g 


m 


vap,g 


Y 


m  +  m  .  +  Am  m  +  m  . 

vap,g  air,g  vap,g  vap,g 


vap,g 


At 

AmvaP'g  mairg  ^ 

At  mg(mg  +  Am  vapg)  mg 


hi 


r,. 


or 


(41) 


(42) 


(43) 


Pg  , 


oc  o  Y  —  ex  — —fh  Y 

LA'g*  g±vaP,g  ** g  ^  rnvap,g±  air ,g 

mg 


m. 


a„ 


V„ 


W) 


111 


(i  -y  ) 

y  \  vap,g ) 

=  M  (l  -  Y  ) 

vap,g  \  vap.g  / 

and  equation  (38)  becomes 


(44) 


dt 


(agPgYvaP,g )  +  V  {otgpgy,Yvaihg )  =  v.(ag  r  V  YvaP'g )  +  MvaP'g  (l  -  YvaP'g ) 


(45) 


In  the  above  equation  M  was  assumed  to  be  known.  An  expression  for 
M  (m  =  - Md),  the  volumetric  mass  rate  of  droplet  evaporation,  will  be 
derived  later. 


Droplet  Balance  Equation 

Droplet  evaporation  is  modeled  by  means  of  the  Unifonn  Temperature  Model  [22,23]. 
This  model  is  based  on  the  assumption  of  a  homogeneous  internal  temperature 
distribution  in  the  droplet  and  phase  equilibrium  conditions  at  the  surface.  The 
analytical  derivation  of  this  model  does  not  consider  contributions  to  heat  and  mass 
transport  through  forced  convection  by  the  gas  flow  around  the  droplet.  Since 
diffusive  time  scales  in  the  surrounding  gas  phase  are  much  smaller  than  that  in  the 
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droplet  fluid,  a  quasi-stationary  description  of  the  gas  phase  is  applied.  The  general 
assumptions  for  the  model  are: 

1 .  The  droplets  are  basically  spherical 

2.  The  exchange  of  mass,  momentum  and  energy  are  quasi-steady 

3.  The  transfer  coefficients  of  mass,  momentum  and  energy  are  expressed 
empirically,  independent  of  the  neighboring  droplets. 


Droplet  Evaporation  Sub-Model 


Due  to  the  difference  in  temperature  and  concentration  between  the  surface  of  a  liquid 
droplet  and  its  surrounding  gas,  a  droplet  may  eventually  transform  into  the  gas-phase. 
This  process  is  called  droplet  evaporation.  The  evaporation  of  droplets  is  essentially  a 
complex  heat  and  mass  transfer  process.  Heat  transfer  first  raises  the  droplet’s 
temperature  until  it  reaches  the  wet  bulb  temperature.  Then,  all  the  heat  transferred  to 
the  liquid  is  used  in  evaporation  and  the  temperature  of  the  rest  of  the  liquid  phase 
remains  unchanged,  this  is  what  is  known  as  the  classical  quasi-steady  uniform 
temperature  droplet  model.  The  evaporation  rate  is  commonly  expressed  as 


dmd 

dt 


(46) 


where  mci  is  the  mass  of  the  liquid  droplet  and,  md  ,  the  mass  flux,  corrected  using  the 
Frossling  correlation  and  based  on  the  classical  droplet  vaporization  model  [11,24], 
Using  reference  values  for  variables  fluid  properties  based  on  the  1/3  rule  of  Sparrow 
and  Gregg  [25],  an  integration  of  the  radially  symmetric  differential  equations  yields 
and  expression  for  the  transport  fluxes.  The  mass  flux  is  given  by: 


md  =  cfm  md  =  cfm 


2nDp,Tg  in 


1-7 


1-7 


vap,s  J 


(47) 


in  this  equation  the  A  indicate  that  the  variable  is  evaluated  at  the  reference 
temperature  and  mass  fraction,  defined  as 


T=—T  +—T 

vap,g  ^  vap,s 

Y  =  — 7  +  —7 

vap,g  ^  vap,s 


(48) 


(49) 


where  T  and  7  are  the  temperature  and  mass  fraction  of  the  vapor  at  the  surface 
of  the  droplet.  Thus  pg  and  F,  are  the  density  and  diffusivity  of  the  gas  evaluated  at 
the  reference  temperature  and  mass  fraction. 


In  the  Unifonn  Temperature  Model,  the  temperature  at  the  droplet  surface  is  basically 
equal  to  that  of  the  droplet: 


Tvap,s  =  Td  (50) 

On  the  other  hand,  the  vapor  concentration  on  the  surface  of  the  droplet  is  found  using 
the  exponential  law  of  Cox-Antoine  [22],  which  states  that: 
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X  ( P  1  =  r sat 

1  vap,s  sat )  p 

with  the  saturation  pressure  Psat  for  a  droplet  at  temperature  Td  obtained  from: 


psa,(Td)  =  e  '  J  (52) 

where  A,  B  and  C  are  specific  values  for  the  droplet  liquid  under  consideration.  Thus 


Y  = 

vap,s 


X  MW 

vap,s  vapor 


vap’s  X  MW  +  (\-  X  \MW . 

vap,s  vapor  vap,sj^  air 

The  liquid  droplet  takes  its  energy  from  the  gas,  which  increases  the  liquid 
temperature.  Once  enough  energy  has  been  supplied  to  overcome  the  latent  heat  of 
evaporation,  evaporation  is  initiated.  This  can  be  mathematically  expressed  as: 

Kap,g,s  ~  K,s  =  AKapor  (Td,S  )  (54) 

The  equation  for  the  heat  transfer  to  the  liquid  droplet  is  given  by: 

dicT.)  . 

™d  +  Qconv,  (55) 

where  hd  is  the  enthalpy  of  the  liquid  droplet,  and  Qmnv,s  and  Qevap  s  are  the  convection 
and  evaporation  heat  transfer  rates,  respectively,  written  as: 

Q<~,,=rt>2f(T„-Td)  (56) 

where  (3 *  in  equation  (56)  is  the  corrected  convective  heat  transfer  coefficient  defined 


p  =cfh/s = cjhj  - :  tai 


'KMtJ  (58) 

The  correction  factors  cfm  and  cjh  [11]  account  for  convective  mass  and  heat  transport 
and  are  computed  from: 

cfm  =  1  +  0.276R4/2  Scf  (59) 


cfm  =  1  +  0.276ReJ  Scf  (59) 

cfh  =  1  +  0.276Rerf1/2Pr^1/3  (60) 

where  Red,  Scd,  and  Prd  are  the  Reynolds,  Schmidt,  and  Prandtl  numbers  respectively, 
defined  as: 


Red  = 


P*  va-vd  Dd 


(62) 


Supersonic  Turbulent  Mixing  and  Evaporation 


19 


Pr,= 


LL  C 

_  P,g 


(63) 


X  and  f  are  the  conductivity  and  diffusivity  of  the  water  vapor  evaluated  at  the 

reference  temperature.  The  convection  tenn  in  equation  (55)  represents  the  amount  of 
energy  due  to  the  interaction  between  the  droplet  and  gas  phases,  and  is  thus  equal  in 
magnitude  to  the  Sh,g  tenn  in  the  energy  equation  of  the  gas  phase,  that  is: 


Sh.g  ~  Qconv,s  (64) 

Note  that  the  energy  equation  for  the  droplet  can  also  be  written  as 
ditnJh)  _  d(hd)  ,  d{md) 

dt  dt  dt  (65) 

= '»  AK{t, J+  kD2/3(t,.,  ~Tt)+  mdh „ 

in  this  case  mass  balance  is  accounted  for  in  the  energy  balance. 


Droplet  Diameter  Equation 


A  relation  between  the  mass  and  the  droplet  diameter  can  be  derived  as  follows. 
Starting  with  the  mass  of  a  droplet  of  volume  Vd,  we  have: 

md  =  PtTd=  Pd  ~7~  (66) 

6 

As  the  droplet  vaporizes,  the  rate  of  mass  decrease  is: 


dmd  _  7 iD2  dD 

dt  2  dt 


(67) 


Finally  before  deriving  the  Eulerian  fonnulation  of  the  droplet  mode,  we  define  for 
later  use  the  volumetric  mass  rate  of  droplet  evaporation  as 


M 


d 


md  _  6a d 

/  \  -  "5 

C Vd!ad )  nDi 


(68) 


Droplet  Discrete  Liquid  Phase 

The  droplet  liquid  phase  governing  equations  are  derived  by  re-formulating  the  droplet 
Lagrangian  model  into  the  Eulerian  framework.  This  is  performed  through  the  use  of 
the  material  derivative. 


Mass  Conservation 


In  the  Eulerian  approach  an  additional  continuity  equation  is  needed  for  the  droplet 
phase.  For  the  mass  conservation  we  have 


a_ 

dt 


(txdpd)+V.{ad\ipd)  = 


V. 


Pturb,d 


\ 

Vad 


+  Md 


V  S Cturb,d 


J 


(69) 
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Again  we  have  a  tenn  on  the  RHS,  which  is  similar  in  fonn  to  diffusion  but  which 
really  accounts  for  the  dispersion  of  the  droplet  in  the  gas. 

Momentum  Conservation 

In  the  Lagrangian  approach,  a  momentum  equation  for  every  droplet  is  derived 
through  a  force  balance  and  is  given  as: 


D(pdVdya) 


-Vd'VP  +  P1Vdi  +  FD  +  mdvl 


In  the  Eulerian  approach  the  droplets  are  no  longer  tracked  individually,  but  the 
ensemble  average  is  used  i.e.  using  the  droplet  (liquid)  volume  fraction: 


When  the  Lagrangian  equation  is  summed  over  all  droplets  and  the  volume  fraction 
substituted,  the  Eulerian  momentum  equation  for  the  droplet  phase  is  found  to  be: 


(^P^d)  +  V.(a,vdp,)  =  -^VJP  +  FB  +Fd  +  M  d\A 


dt  a 
where 


Fd  =-Fd 

rd  rg 


F,?=«/A/g  (74) 

For  turbulent  flow,  an  extra  term  is  added  on  the  right  hand  side  of  Eq.  (74),  which 
should  not  be  interpreted  as  a  diffusion  tenn  since  the  droplets  can  be  regarded  as 
“solid  particles”  and  thus  non-diffusive  entities.  Instead,  this  term  represents  the 
dispersion  of  the  droplets  in  the  gas  due  to  turbulent  interaction,  and  this  dispersion 
tenn  have  a  similar  fonn  as  the  diffusion  term  in  the  momentum  equation  of  the  gas 
phase,  thus  the  momentum  equation  for  the  droplet  phase  becomes: 

jt  (adPd  vd ) + V.(«,  v<lp(1)  =  -adVP  +  ^ 

V-k/w  [Vvd+Vvdr])+  otdpdg  +  FD  +Md\d 

Droplet  Diameter  Equation 

A  relation  between  M  d  and  m  d  can  be  readily  derived,  defining  the  volumetric 
mass  rate  as  [22,26], 
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Thus 
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dt  3otdpd 

In  order  to  derive  the  transport  equation  for  the  droplet  diameter  we  first  write 


V {a dydpdD)  =  DV.(adydpd)+  ad\Apd  ■  VD 

noting  that 


dD 

dt 


—  +  vdVD 

dt  d 


we  get: 


(77) 
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(79) 
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To  account  for  turbulence  the  above  equation  is  written  as: 


(80) 
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Energy  Conservation 

A  transport  equation  for  the  droplet  temperature  can  also  be  derived  via  and  energy 
balance  in  the  same  manner  as  the  momentum  equation,  this  equation  is  written  as: 

^(A/A/0+v-(A/AAA) 

=  K  |:k/A/)+/f/V-k/A/Vd)+ A/A/  I^J+A/AAd-V^J  (82) 


:  A/ A, 
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dt  p'd  dt 


A+V.|  Atm  Va, 

i  ^Cturb,d 


Following  the  derivation  of  Faeth  [22]  and  Wittig  [26]  in  relation  to  droplet  heating 
and  evaporation  using  the  so  called  “Uniform  Temperature  Model”,  the  balance 
equations  for  the  droplet  mass  and  temperature  are 
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„  ,e,  v  1  Z.JUJ/1  refc  m  1 

e  *<•«  - 1  e  g  - 1 

In  the  above  equation  the  first  tenn  on  the  RHS  accounts  for  heat  exchange  with  the 
surrounding  and  the  second  term  on  the  RHS  accounts  for  the  evaporation  latent  heat. 

The  energy  equation  is  now  rewritten  as: 
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To  account  for  turbulence  the  above  equation  is  rewritten  as: 
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(87) 


nD2[f(Tg-Td)+Md(Ahv  +  hd) 

where  pturhd  is  the  turbulent  viscosity  in  the  disperse  phase  computed  using  the 
algebraic  model  described  next. 


Turbulence  Model  in  the  Secondary  Dispersed  phase 


Many  workers  [3,13,14]  have  successfully  used  the  closure  hypotheses  for  the 
turbulent  mass,  momentum,  and  heat  transfer  fluxes  used  in  this  work.  Moreover,  the 
turbulent  viscosity  of  the  disperse  (droplet)  phase,  ptUrb,d,  is  modeled  using  the 
approach  of  Melville  and  Bray  [27]: 
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A  fundamental  assumption  of  this  approach  is  the  dependence  of  the  turbulent 
viscosity  of  the  droplet  phase  on  local  mean  flow  properties.  As  for  the  ratio  of  the 
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turbulent  kinetic  energies  of  the  dispersed  and  gas  phase  it  is  calculated  following  the 
approach  of  Kramer  [28]  as: 


kd  _  1 

kg  1  +  co2t2 


(90) 


where  t  is  the  particle  relaxation  time,  and  co  is  the  frequency  of  the  particle  response 
which  are  given  by: 


z_lPdD2  1 

18  pg  v  1  +  0.133  Re  °'687 

with  the  following  characteristic  macroscopic  length  scale  for  turbulence 


(91) 


L 


x 


(92) 


Since  in  general  the  droplets  do  not  follow  the  motion  of  the  surrounding  fluid  from 
one  point  to  another  it  is  expected  for  the  ratio  kti/kg  to  be  different  from  unity  and  to 
vary  with  particle  relaxation  time  t  and  local  turbulence  quantities.  For  the  turbulent 
Schmidt  number  of  the  droplet  phase,  Sclurbd,  Kramer  [28]  suggests  a  value  of  0.3. 
However  in  more  recent  work  [4]  it  was  found  to  be  particle  size  dependent  and  a 
value  of  0.7  is  used  in  this  work  ( Scturbd  =  0.7).  For  the  turbulent  Prandtl  number  a 


value  of  0.85  was  chosen  (Pr,„,.M  =  0.85  ). 


Numerical  Approach 

A  review  of  the  above  differential  equations  reveals  that  they  are  similar  in  structure. 
If  a  typical  representative  variable  associated  with  phase  (k)  is  denoted  by  <])(k),  the 
general  fluidic  differential  equation  may  be  written  as: 

d(cc(k)o(k)rtk)} 

-A — kJW+V.(/ytVt)f)=  V.(a(AT,i'1V^))+  dk)Qik)  (93) 

where  the  expression  for  r<k)  and  Q(k)  can  be  deduced  from  the  parent  equations. 


Discretization  Procedure 

In  the  previous  sections  the  differential  equations  governing  multi-fluid  flow 
phenomena  of  spray  mixing  and  evaporation  were  outlined.  The  task  now  is  to  present 
the  Finite  Volume-based  numerical  solution  algorithm  for  solving  these  equations. 
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Discretization  of  the  general  fluidic  conservation  equation 


The  general  conservation  equation  (93)  is  integrated  over  a  finite  volume  to  yield: 

rr 

^q+  j] 

*  o  (94) 

=  Jjv.(eraTaVfA,)/Q  +  JjV^A/Q 

n  n 

Where  Q  is  the  volume  of  the  control  cell.  Using  the  divergence  theorem  to  transform 
the  volume  integral  into  a  surface  integral  and  then  replacing  the  surface  integral  by  a 
summation  of  the  fluxes  over  the  sides  of  the  control  volume,  equation  (91)  is 
transformed  to: 


d(a(k)p(k>(/)(k,n)  mr.  ...  ... 

v  — J-+  Z(J  T  +  oc(k)Q,k)Q  (95) 

'y'  nb=e,w,n,s,t,b 

where  Jf  D  and  Jf  c  are  the  diffusive  and  convective  fluxes,  respectively.  The 
discretized  form  of  the  diffusive  flux  along  an  east  face  is  given  by: 


r<*>0 


.(i-)-p(i-) 


=  -a  T 


(4>\ 


'T 


(96) 


where  the  over  bar  denotes  a  value  obtained  by  interpolation,  Ke  is  a  space  vector 
defined  as, 


k,  =  (A.,  -  (?d) y,  =  <i  +  <j  +  r;k  (97) 

and  y  is  a  scaling  factor  given  by  [29], 


7e  = 


Sede 

Sede 


(98) 


and  nfand  df  are  the  contravariant  and  covariant  unit  vectors  respectively.  The 
discretized  convective  flux  along  an  east  face  is  given  by: 


fk)C  =  (efV,uW  ■  S)  =  (c/k)p(-k)U(k))e  <j>f  =  Of  (99) 

where  Se  and  Ce  are  the  surface  vector  and  convection  flux  coefficient  at  cell  face  'e', 
respectively.  As  can  be  seen  from  Eq.  (99),  the  accuracy  of  the  control  volume 
solution  for  the  convective  scalar  flux  depends  on  the  proper  estimation  of  the  face 
value  <t>e  as  a  function  of  the  neighboring  0  nodes  values.  Using  some  assumed 
interpolation  profile,  4>e  can  be  explicitly  fonnulated  in  tenns  of  the  nodal  values  by  a 
functional  relationship  of  the  fonn: 

A  =  /fc)  (100) 

where  c|)nb  denotes  the  <|)  values  at  the  neighboring  nodes.  The  interpolation  profile  may 
be  one-dimensional  or  multi-dimensional  of  low  or  high  order  of  accuracy.  The  higher 
the  order  of  the  profile  is,  the  lower  numerical  diffusion  will  be.  However,  the  order  of 
the  profile  and  its  dimensionality  do  not  eliminate  numerical  dispersion.  To  minimize 
this  error,  limiters  on  the  convective  flux  should  be  imposed  to  enforce  monotonicity 
[30].  For  more  details  the  reader  is  referred  to  [31]. 
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After  substituting  the  face  values  by  their  functional  relationship  relating  to  the  node 
values  of  (j),  Eq.  (95)  is  transformed  after  some  algebraic  manipulations  into  the 
following  discretized  equation: 

=£<«+*?>  (ioi) 

NB 

where  the  coefficients  Ap]  and  depend  on  the  selected  scheme  and  Bp]  is  the 
source  term  of  the  discretized  equation  .  In  compact  form,  the  above  equation  can  be 
written  as 

y  Aa)dtk)  +  B[k) 

Z-i^NBYNB  +  DP 

^  =  - - -  (102) 


Discretization  of  the  fluidic  momentum  equation 


The  discretization  procedure  for  the  momentum  equation  yields  an  algebraic  equation 
of  the  fonn: 


A<*>u<‘> 


=  2><> 


D1 

NB  "r  DP 


(*) 


■otpapj  r(p)+arYdgw(f‘P  ■ 


u 


NB(P ) 


m^k 


(103) 


In  the  above  equation,  the  inter-phase  term  is  written  out  explicitly  to  show  the  strong 
coupling  among  the  momentum  equations  of  the  different  fluids.  This  is  different  from 
the  spatial  coupling  that  exists  among  the  neighboring  velocities  of  the  same  fluid.  To 
improve  the  overall  convergence  and  robustness  of  the  algorithm,  the  inter-phase 
source  term  can  be  linearized  and  the  linear  part  treated  implicitly  so  that: 
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Where  now 

AW^A(;,  +  QpX(gpm))  (105) 

m^k 

For  later  reference,  the  discretized  fonn  of  the  momentum  equation  is  written  in  the 
following  two  fonns: 

Up  >  =  HP p  [uw]-  Nk,D<k)Vp(P)  (106) 

where  the  body  force  and  inter-phase  tenns  are  absorbed  in  the  Bj,k) source  term  within 
the  HPp  [u(k)  ]  term,  or  as 

Up  >  =HIp  [u(i)]+D(/)X(g(fe”)<m))  (107) 

m^k 

where  the  body  force  and  pressure  gradient  terms  are  absorbed  in  the  Bpk  1  source  tenn 
within  the  HIP  [u<kl]  term. 
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Discretization  of  the  fluidic  mass  conservation  equation 

The  fluidic  mass-conservation  equation  can  be  viewed  as  a  fluidic  volume  fraction 
equation: 

dp  =  HP  [dk)\  (108) 

or  as  a  fluidic  continuity  equation  to  be  used  in  deriving  the  pressure  correction 
equation: 

(nt(k)  dkA-(fy(k)  dkA°ld  r  l 

P  PP  &  ‘  P''  Q  +  )PWuW-s]  = M(  )  (109) 

where  the  A  operator  represents  the  following  operation: 

A  p[0]=  E©/  (HO) 

f=NB(P ) 

Discretization  of  the  fluidic  energy  equation 

The  discretization  of  the  energy  equation  follows  that  of  the  general  fluidic  scalar 
equation.  The  only  difference  is  the  one  pertaining  to  the  discretization  of  the 
additional  source  terms.  Since  a  control  volume  approach  is  followed,  the  integral  of 
these  sources  over  the  control  volume  appears  in  the  discretized  equation.  By  using  the 
divergence  theorem,  the  volume  integral  is  transformed  into  a  surface  integral  and  the 
resultant  discretized  expressions  evaluated  explicitly. 


Solution  Procedure 


The  number  of  equations  describing  an  n-fluid  flow  situation  is:  n  phasic  momentum 
equations,  n  phasic  volume  fraction  (or  mass  conservation)  equations,  a  geometric 
conservation  equation,  and  for  the  case  of  compressible  flow  additional  n  auxiliary 
pressure-density  relations.  Moreover,  the  variables  involved  are  the  n  phasic  velocity 
vectors,  the  n  phasic  volume  fractions,  the  pressure  field,  and  for  a  compressible  flow 
an  additional  n  unknown  phasic  density  fields.  In  the  current  work,  the  n  momentum 
equations  are  used  to  calculate  the  n  velocity  fields,  n-1  volume  fraction  (mass 
conservation)  equations  are  used  to  calculate  n-1  volume  fraction  fields,  and  the  last 
volume  fraction  field  calculated  using  the  geometric  conservation  equation 

dn)  =  \~Y,oc(k)  (111) 

k^n 

The  remaining  volume  fraction  equation  can  be  used  to  calculate  the  pressure  field  that 
is  shared  by  all  phases.  However,  instead  of  using  this  last  volume  fraction  equation, 
in  the  class  of  Mass  Conservation  Based  Algorithms  (MCBA)  the  global  conservation 
equation  is  employed,  i.e.  the  sum  of  the  individual  mass  conservation  equations,  to 
derive  a  pressure  correction  equation  as  outlined  next. 

The  Pressure  Correction  Equation 

To  derive  the  pressure-correction  equation,  the  mass  conservation  equations  of  the 
various  phases  are  added  to  yield  the  global  mass  conservation  equation  given  by: 


Supersonic  Turbulent  Mixing  and  Evaporation 


28 
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Q  +  Ap(awp^uw.s) 


V  M[k)  =0 


(112) 


In  the  predictor  stage  a  guessed  or  an  estimated  pressure  field  from  the  previous 
iteration,  denoted  by  P°,  is  substituted  into  the  momentum  equations.  The  resulting 
velocity  fields  denoted  by  u ' k  1  which  now  satisfy  the  momentum  equations  will  not, 
in  general,  satisfy  the  mass  conservation  equations.  Thus,  corrections  are  needed  in 
order  to  yield  velocity  and  pressure  fields  that  satisfy  both  equations.  Denoting  the 
corrections  for  pressure,  velocity,  and  density  by  P',  u<k)  ,  and  p<k)  respectively,  the 
corrected  fields  are  written  as: 


P  =  P°  +  P',u(k)  =  +  u  W’,p(k)  =  pa)°  +  p(k1  (113) 

Hence  the  equations  solved  in  the  predictor  stage  are: 

u(;}*  =  Hp[uw]-^)°D?)VpP°  (114) 

While  the  final  solutions  satisfy 

u«  =  Hp  [ua)]  -  a(k  D  fVpP  (115) 

Subtracting  the  two  equation  sets  ((115)  and  (114))  from  each  other  yields  the 
following  equation  involving  the  correction  terms: 

<»  =  Hp[ua ’]  -  a[k)°D[k)V  PP'  (116) 


Moreover,  the  new  density  and  velocity  fields,  p<k)  and  u(k),  will  satisfy  the  overall 
mass  conservation  equation  if: 
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Linearizing  the  (p(klu(k))  term,  one  gets 
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(p  +p  J(u  +u  j=p  u  +p  u  +p  u  +p  u  (118) 

Substituting  equations  (118)  and  (116)  into  equation  (117),  rearranging,  and  replacing 
density  correction  by  pressure  correction,  the  final  form  of  the  pressure-correction 
equation  is  written  as: 
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(119) 


The  corrections  are  then  applied  to  the  velocity,  pressure,  and  density  fields  using  the 
following  equations: 

=u(k)°-a(k)°D(k)VpP'  P*=P°  +  P'  p(kf  =  pik)°  +  C(k)P'  (120) 

Numerical  experiments  using  the  above  approach  to  simulate  air-water  flows  have 
shown  poor  conservation  of  the  lighter  fluid.  This  problem  is  considerably  alleviated 
by  normalizing  the  individual  continuity  equations,  and  hence  the  global  mass 
conservation  equation,  by  means  of  a  weighting  factor  such  as  a  reference  density  p(k) 
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(which  is  fluid  dependent).  This  approach  has  been  adopted  in  solving  all  problems 
presented  in  this  work  (see  [32]  for  details). 

The  MCBA-SIMPLE  Algorithm 

The  overall  solution  procedure  is  an  extension  of  the  single-phase  SIMPLE  algorithm 
into  multi-phase  flows.  Since  the  pressure  correction  equation  is  derived  from  overall 
mass  conservation,  it  is  denoted  by  MCBA-SIMPLE  [31].  The  sequence  of  events  in 
the  MCBA-SIMPLE  is  as  follows: 

1 .  Solve  the  fluidic  momentum  equations  for  velocities. 

2.  Solve  the  pressure  correction  equation  based  on  global  mass  conservation. 

3.  Correct  velocities,  densities,  and  pressure. 

4.  Solve  the  fluidic  mass  conservation  equations  for  volume  fractions. 

5.  Solve  the  fluidic  scalar  equations  (k,  8,  T,  Y,  D,  etc. . .). 

6.  Return  to  the  first  step  and  repeat  until  convergence. 
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Physical  Properties 

In  the  test  problems  that  follow  the  material  properties  used  for  the  water  droplets, 
water  vapor  and  air  are  as  given  by  [33]: 

MWair  =  28.97 kg/  kg  mole 

juair  =  6. 109x10  6  +  4.604x1 0~8  7  - 1 .05 1x10 ~UT2  Kg  /  m.s 
Xair  =  3.227x1 0“3  +  8.3894x1 0~57-  1.9858x1 0~8T2J/m.s.K 
Prmr  =  0.8 1 5  -  4.958x1 0  4 7  +  4.5 14x1  0~7  72  forT  <  600 K 
Pr,!r  =  0.647  +  5.5x10 ~ST  forT  >  600K 
MWwater  =18.015 kg  I  kg  mole 

cp  v  =8137-  37.34x7  +  0.07482x72  -  4.956x10  5T2J / kg.K  (121) 

juv  =4.07x1  O'8 T  -3.077x1 0^6  Kg  I  m.s 

Av  =  1 .024x1 0“2  -  8.2  lxl  0“6 T  + 1 .4 lxl 0”7 T2  -  4.5  lxl  0“u T3  J  /  m.s.K 
Ahv  =2.257x1 06  +  2.595x1 02  {373.15 -T)j  /  kg 
pL  =  997 kg  /  m2 
cL  =  41 84  J  /  kg.K 
X,  =  0.653  \J/m.s.K 

It  is  further  assumed  that  the  Lewis  number  is  equal  to  unity  and  the  diffusivity  (r)  is 
calculated  with  the  density  evaluated  at  the  reference  temperature. 

The  use  of  the  above  relations  and  materials  are  for  convenience  and  should  not  be 
taken  as  restrictions  related  to  the  code. 
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Results  and  Discussion 


In  what  follows,  solutions  to  three  two-dimensional  two-phase  flow  problems  are 
presented  and  discussed.  The  physical  situations  for  these  problems,  depicted  in  Figs. 
4(a)-4(c),  represent  a  rectangular  duct  in  which  air  enters  with  a  uniform  free  stream 
velocity  U,  while  fuel  (water  is  used  here)  mixed  with  air  is  injected  through  a  nozzle 
2mm  in  diameter  either  in  the  stream- wise  direction  (Fig.  1(a)),  or  in  the  cross  stream 
direction  (Figs.  1(b)  and  1(c)).  The  length  of  the  domain  is  L  (L=l  m)  and  its  width  is 
W  (W=0.25m).  Turbulent  flow  results  are  generated  using  the  k-e,  k-co,  and  SST 
turbulence  models.  Moreover,  the  required  temperature  dependent  properties  of  the 
carrier  gas  (air)  and  vapor  species  are  as  given  above. 


Injection  in  the  stream-wise  direction 


The  physical  domain  depicted  in  Fig.  4(a)  is  subdivided  into  120x89  non-unifonn 
control  volumes.  The  fuel  is  injected  through  9  uniform  control  volumes  (each  of 
width  0.002/9  m)  with  equal  x-component  of  velocity  (=20%  of  the  free  stream 
velocity)  and  different  injection  angles  (varying  unifonnly  from  -45°  to  45°  as  shown 
in  Fig.  4(a)). 

To  validate  the  implemented  evaporation  model,  the  two-phase  solution  is  compared 
against  the  single  phase  solution  for  the  turbulent  supersonic  case  (Miniet=2)  using  the 
k-8  model.  The  single  phase  problem  corresponds  to  exactly  the  same  two-phase 
problem  but  with  the  injected  water-air  mixture  replaced  by  air.  Results  generated  are 
displayed  in  Fig.  5.  In  these  plots,  cross-stream  profiles  for  the  x-component  of  the  gas 
velocity  (Fig.  5(a)),  the  gas  temperature  (Fig.  5(b)),  the  gas  density  (Fig.  5(c)),  and  the 
gas  turbulent  kinetic  energy  (Fig.  5(d))  at  three  axial  stations  are  compared.  As 
depicted  in  Fig.  5(a),  the  two-phase  U-velocity  profiles  indicate  lower  values  than  the 
single  phase  profiles,  which  is  physically  expected  as  the  presence  of  particles  injected 
at  lower  velocity  than  the  gas  phase  creates  drag  that  slows  down  the  gas  phase.  In 
Fig.  5(b),  the  lower  temperature  values  of  the  two-phase  flow  are  also  in  accord  with 
what  is  expected  since  the  droplet  phase  absorbs  heat  as  it  evaporates  and  thereby 
decreases  the  temperature  of  the  gas  phase.  Moreover,  the  higher  density  in  the 
presence  of  droplets  (Fig.  5(c))  is  due  to  the  water  evaporating  into  air.  As  shown,  the 
difference  increases  with  distance  as  the  mass  fraction  of  water  vapor  increases.  The 
turbulent  kinetic  energy  profiles  displayed  in  Fig.  5(d)  indicate  higher  turbulence  level 
in  the  presence  of  droplets,  which  is  physically  correct. 

The  size  of  the  grid  used  in  the  computations  to  obtain  grid  independent  solutions  was 
selected  after  comparing  results  generated  on  grids  of  increasing  density.  As 
mentioned  above,  a  final  grid  of  size  121x91  grid  points  was  chosen  as  results 
generated  with  this  grid  are  almost  identical  to  the  ones  generated  using  a  grid  of  size 
162*  111  grid  points.  This  is  clearly  seen  in  Fig.  6  where  cross-stream  profiles  of  the  x- 
component  of  the  gas  velocity  (Fig.  6(a)),  gas  temperature  (Fig.  6(b)),  vapor  mass 
fraction  in  the  gas  phase  (Fig.  6(c)),  and  droplet  temperature  (Fig.  6(d))  at  three  axial 
stations  are  compared. 
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Having  established  physical  and  numerical  credibility,  the  problem  is  solved,  as 
described  next,  for  subsonic  and  supersonic  conditions. 


Subsonic  evaporation  and  mixing  of  water  droplets  in  air 


For  the  physical  situation  depicted  in  Fig.  4(a),  the  Mach  number  and  temperature  of 
the  air  at  inlet  to  the  domain  are  taken  to  be  0.2  (Mair,iniet=0.2)  and  900  K,  respectively. 
The  mixture  of  air  and  droplet  are  injected  into  the  domain  with  a  temperature  of  300 
K  with  the  volume  fraction  of  the  water  in  the  injected  air- water  mixture  being  10'“. 
The  velocity  of  the  injected  mixture  is  given  as  (there  is  no  particular  reason  for  the 
choice  of  this  velocity  profile): 


d,  injected 


0.2*  v  .  ,  i  +  0.2*  v  .  ,  sin 


II  g,  inlet 


II  g,  inlet 


n  k-\n 


=  0.3 *  v  ..  i  +  0.3*  v  . ,  sin 


g,  injected 


\\g, inlet 


\\g, inlet 


A  4  4 
n  k  - 1  n 


4  4  4 


j  1  <  k  <  9 


j  1  <  k  <  9 


(122) 


With  this  injection  velocity  profile,  there  is  a  slip  between  the  two  injected  phases  and 
the  magnitude  of  the  injected  velocity  increases  away  from  the  center  of  the  jet. 
Moreover,  with  this  velocity  profile  and  volume  fraction  a  total  of  0.479Kg/s  of  fuel 
are  injected  into  the  domain,  which  corresponds  to  an  overall  fuel-air  ratio  of  0.04124 
kgfaei/kgair-  This  fuel/air  ratio  is  at  the  high  side  of  typical  values  used  in  gas  turbine 
applications. 


The  problem  is  first  solved  assuming  the  flow  to  be  laminar  and  then  turbulent  using 
the  k-e  model  and  results  are  displayed  in  Figs.  7-9.  In  Fig.  7,  contour  plots  for  the 
laminar  case  are  presented.  In  Figs.  7(a)  and  7  (b)  the  droplet  velocities  and  volume 
fraction  fields  are  depicted.  As  shown,  due  to  the  higher  air  velocity,  the  spreading  of 
injected  water  is  low  and  droplets  quickly  align  with  the  air  velocity.  As  expected,  the 
gas  temperature  (Fig.  7(c))  decreases  along  the  centerline  of  the  domain  because  of  the 
droplet  evaporation.  This  decrease  in  temperature  is  associated  with  an  increase  in 
density.  The  distribution  of  water  vapor  in  the  gas  phase  is  displayed  in  Fig.  7(d)  while 
the  pressure  field  is  depicted  in  Fig.  7(e).  The  mass  fraction  of  the  water  vapor  in  the 
gas  phase  increases  as  the  mixture  moves  downstream  in  the  channel  due  to  the 
increase  in  the  evaporated  amount  with  distance,  which  is  physically  plausible. 
Computations  have  shown  that  26.953%  of  the  water  droplets  entering  the  domain  are 
transformed  into  water  vapor. 


Turbulent  flow  results  for  the  same  configuration  and  physical  properties  using  the  k-8 
turbulence  model  are  displayed  in  Fig.  8.  Due  to  the  higher  mixing  in  turbulent  flows, 
the  droplets  align  faster  with  the  gas  field  resulting  in  less  spreading  of  the  droplets  in 
comparison  with  the  laminar  case,  as  depicted  by  the  droplet  velocity  and  volume 
fraction  fields  presented  in  Figs.  8(a)  and  8(b).  For  the  same  reason,  a  higher  decrease 
in  temperature  is  obtained  along  the  central  portion  of  the  domain  (Fig.  8(c))  where  the 
water  vapor  is  concentrated  (Fig.  8(d)).  Unlike  the  laminar  case,  the  pressure  decreases 
in  the  stream-wise  direction  (Fig.  8(e)).  To  be  mentioned  is  the  fact  that  26.947%  of 
the  water  droplets  entering  the  domain  are  transformed  into  water  vapor.  This  is  0.06% 
lower  than  the  value  obtained  in  the  laminar  case  and  can  be  explained  as  follows:  The 
two  factors  controlling  the  rate  of  evaporation  are  the  residence  time  and  rate  of 
mixing.  As  the  rate  of  increase  in  the  x-component  of  the  droplet  velocity  is  lower  for 
laminar  flow  than  for  turbulent  flow  (due  to  the  lower  mixing  rate),  the  residence  time 
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is  higher  for  laminar  than  for  turbulent  flows.  Consequently,  in  the  laminar  case, 
evaporation  occurs  over  a  longer  period.  However,  this  increase  in  residence  time  is 
associated  with  lower  mixing  and  heat  transfer  rates,  which  is  exactly  opposite  to  what 
takes  place  in  the  turbulent  case.  These  two  competing  and  opposite  effects  seem  to 
balance  each  other  and  tend  to  keep  the  rate  of  evaporation  nearly  the  same  for  laminar 
and  turbulent  flows,  even  though  the  profiles  of  the  various  fields  are  different. 

The  decrease  in  the  droplet  diameter  and  increase  in  droplet  temperature  along  the 
centerline  of  the  channel  is  depicted  in  Fig.  9  for  the  turbulent  flow  case.  As  shown, 
the  rate  of  increase  in  temperature  is  high  in  the  early  part  of  the  channel  where  90% 
of  the  total  increase  occurs  over  almost  25%  of  the  channel  length  due  to  the  relatively 
low  gas  velocity.  This  increase  in  temperature  is  associated  with  a  local  decrease  in  the 
rate  of  evaporation,  which  explains  the  slight  decrease  in  the  rate  at  which  the 
diameter  decreases. 


Supersonic  evaporation  and  mixing  of  water  droplets  in  air 

For  the  same  physical  situation  depicted  in  Fig.  4,  the  air  Mach  number  at  inlet  to  the 
domain  which  is  set  to  2  CMajrjn|et=2)  while  the  values  of  all  other  parameters  are  the 
same  as  in  the  subsonic  case.  Results  are  presented  in  Figs.  10-14.  In  Fig.  10, 
generated  results  for  the  laminar  case  are  displayed.  Similar  results  for  turbulent  flow 
using  the  k-8,  k-co,  and  SST  models  are  depicted  in  Fig.  11,  12,  and  13,  respectively.  In 
all  Figures,  field  plots  of  particle  velocity  (Figs.  11(a), 12(a),  and  13(a)),  particle 
volume  fraction  (Figs.  11(b), 12(b),  and  13(b)),  gas  temperature  (Figs.  11(c), 12(c),  and 
13(c)),  mass  fraction  of  water  vapor  in  the  gas  phase  (Figs.  11(d), 12(d),  and  13(d)), 
and  pressure  (Figs.  1 1(e), 12(e),  and  13(e))  are  presented.  Due  to  their  high  inertia  and 
the  much  lower  viscous  effects  in  the  laminar  flow  case,  the  droplets  are  less  affected 
by  the  gas  and  require  a  longer  distance  to  align  with  the  primary  gas  flow  in 
comparison  with  the  turbulent  flow  cases  (compare  the  particle  velocity  fields  in  Figs. 
10(a)- 13(a)).  In  this  case,  the  two  forces  other  than  inertia  are  the  pressure  gradient 
and  drag.  This  explains  the  completely  different  pressure  field  obtained  in  the  laminar 
case  (Fig.  10(e))  as  compared  to  pressure  fields  obtained  with  the  various  turbulence 
models  (Figs.  1 1(e)-  13(e)).  Moreover,  the  general  trend  in  variation  of  other  variables 
is  similar  and  resembles  the  subsonic  case,  i.e.  the  volume  fraction  of  the  particles 
decreases  in  the  stream-wise  direction,  the  gas  temperature  decreases  along  the  central 
portion  of  the  domain  due  to  the  absorption  of  heat  by  the  droplets  to  evaporate,  and 
the  stream-wise  increase  in  the  mass  fraction  of  the  water  vapor  in  the  gas  phase  as 
more  water  evaporates. 

To  be  noticed  here  are  the  droplet  velocity  profiles  predicted  by  the  various  turbulence 
models  displayed  in  Figs.  ll(a)-13(a).  The  highest  spreading  is  obtained  with  the  k-8 
model  and  the  lowest  spreading  by  the  k-co  model.  The  predictions  of  the  k-co  model, 
unlike  the  k-8  model,  are  well  known  to  depend  on  the  free  stream  values,  ft  seems 
that  the  inlet  turbulence  quantities  used  in  the  simulation  do  not  work  well  with  the  k- 
co  model.  The  inlet  turbulence  intensity  is  taken  to  be  1=5%  and  the  turbulence  length 
scale  L=0.25  (the  height  of  the  domain)  and  k,  8,  and  co  at  inlet  are  computed  as: 
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In  fact,  by  setting  the  value  of  to  at  inlet  to  be  equal  to  that  of  8,  spreading  close  to  the 
one  obtained  in  the  laminar  case  was  predicted.  Therefore  to  be  able  to  reproduce  the 
k-8  predictions,  the  value  of  co  at  inlet  has  to  be  adjusted  through  trial  and  error.  This 
task  has  not  been  undertaken  due  to  time  limitation.  The  SST  model  is  supposed  to 
possess  the  virtues  of  both  the  k-8  model  away  from  walls  and  the  virtues  of  the  k-co 
model  close  to  the  wall  and  by  including  the  cross  diffusion  terms,  it  becomes  free- 
stream  independent.  Moreover,  its  predictions  are  expected  to  lie  between  those  of  the 
k-8  and  k-co  models.  As  expected,  results  presented  in  Fig.  13  possess  these  properties. 
In  fact,  the  highest  turbulent  viscosity  over  the  domain,  based  on  the  above  inlet 
turbulence  quantities,  is  close  to  4  Kg/(m.s)  for  the  k-8  model,  close  to  28  Kg/(m.s)  for 
the  k-co  model,  and  close  to  6  Kg/(m.s)  for  the  SST  model.  Clearly,  the  higher  the 
turbulent  viscosity  is,  the  less  spreading  is  predicted. 

The  increase  in  the  droplet  temperature  and  decrease  in  the  droplet  diameter  along  the 
centerline  of  the  channel,  as  predicted  with  the  k-8  model,  is  depicted  in  Fig.  14.  Due 
to  the  high  speed  of  the  flow,  the  droplet  temperature  does  not  reach  its  limiting  wet 
bulb  value.  Moreover,  the  decrease  in  the  droplet  diameter  is  lower  than  in  the 
subsonic  case.  In  fact,  computations  have  shown  that  only  around  4.9%  of  the  water 
droplets  entering  the  domain  are  transformed  into  water  vapor. 

To  study  the  effect  of  the  various  parameters  on  the  evaporation  rate,  the  SST  model  is 
selected  and  results  are  presented  in  Figs.  15-17.  The  effect  of  the  inlet  gas 
temperature  is  presented  in  Fig.  15.  For  that  purpose,  all  parameters  are  held  constants 
with  the  exception  of  the  inlet  gas  temperature,  which  is  varied  between  300  and  1100 
K.  As  expected,  the  rate  of  evaporation  increases  with  increasing  T„jnicl.  The  highest 
value  achieved  is  around  5.5%,  which  is  still  low  and  clearly  demonstrates  the 
difficulties  one  encounters  in  high  speed  flows.  The  effect  of  the  droplet  temperature 
on  the  rate  of  evaporation  is  displayed  in  Fig.  16.  As  shown,  the  Tdjiniet  is  varied 
between  300  and  370  K,  while  holding  all  other  parameters  constants.  The  highest 
value  used  is  slightly  less  than  the  water  boiling  temperature  at  normal  atmospheric 
pressure.  As  depicted  in  Fig.  16,  the  rate  of  evaporation  increases  with  increasing  the 
droplet  temperature.  It  can  also  be  inferred  from  the  figure  that  it  is  more 
advantageous  to  increase  the  droplet  temperature  than  the  gas  temperature,  as  the 
highest  rate  of  evaporation  achieved  with  an  inlet  droplet  temperature  of  370  k  is  close 
to  14.4%.  Finally,  the  effect  of  the  domain  length  on  the  rate  of  evaporation  is  studied 
and  results  reported  in  Fig.  17.  For  that  purpose,  a  droplet  temperature  of  370  K  is 
selected  and  the  length  of  the  domain  varied  from  0.5  to  2  m  and  computations 
perfonned  holding  all  other  parameters  constants.  As  expected,  the  rate  of  evaporation 
increases  with  increasing  the  domain  length  and  for  the  longest  domain  its  value  is 
close  to  20.5%. 
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Injection  at  Right  Angle  From  One  Nozzle 


The  physical  situation  for  this  problem,  depicted  in  Fig.  4(b),  represents  a  rectangular 
duct  in  which  air  enters  with  a  uniform  free  stream  velocity  U,  while  fuel  (water  is 
used  here)  mixed  with  air  is  injected  from  a  nozzle  on  the  lower  horizontal  wall  of  the 
domain.  The  length  of  the  domain  is  L  (F=l  m)  and  its  width  is  W  (W=0.25m). 
Moreover,  the  required  temperature  dependent  properties  of  the  carrier  gas  (air)  and 
vapor  species  are  given  by  Eq.(121). 

The  Mach  number  and  temperature  of  the  air  at  inlet  to  the  domain  are  taken  to  be  2 
(Mair,iniet=2)  and  900  K,  respectively.  The  mixture  of  air  and  droplet  enters  the  domain 
with  a  temperature  of  300  K  with  the  volume  fraction  of  the  water  in  the  injected  air- 
water  mixture  being  10'2,  the  velocity  of  the  fuel-air  mixture  is  10%  the  gas  velocity  in 
the  main  gas  stream,  with  a  fuel  mass  flow  rate  of  2.398  Kg/s.  An  81x41  grid  system 
is  used  and  the  fuel-air  mixture  is  injected  nonnal  to  the  main  gas  stream  from  a 
nozzle  2mm  in  size  using  five  control  volumes.  Moreover,  results  are  generated  for 
laminar  and  turbulent  flow  using  the  k-e,  k-co,  and  SST  turbulence  models. 

Results  for  this  problem  are  presented  in  Figs.  18-22.  Laminar  flow  plots  are  displayed 
in  Fig.  18  whereas  results  for  turbulent  flow  are  presented  in  Figs.  19-21.  The  higher 
jet  penetration  experienced  by  the  laminar  computations  as  compared  to  the  turbulent 
ones  can  be  explained  as  follows.  As  shown  in  Figs.  1 8(a)-2 1  (a),  the  highest 
turbulence  of  the  gas  phase  is  at  inlet  and  decreases  downstream  due  to  its  interaction 
with  the  nonnal  jet  that  results  in  large  dissipation,  which  is  manifested  by  the  increase 
in  temperature  near  the  nozzle  inlet.  Furthennore,  the  increased  mixing  in  the 
turbulent  simulation  means  that  a  higher  percentage  of  the  droplets  have  evaporated  in 
the  near  nozzle  zone  and  a  lower  momentum  is  available  for  the  droplet  to  penetrate 
the  gas  stream.  These  results  are  consistent  for  all  three  turbulence  models,  with  the  k- 
8  model  yielding  slightly  higher  spreading  than  the  other  two  models.  This  effect  does 
point  out  to  the  importance  that  a  turbulent  gas  phase  can  have  on  the  injected  fuel. 
The  general  trend  of  results  is  also  consistent  with  the  parallel  test  problem  where  the 
effects  of  the  gas  phase  on  the  droplets  being  larger  in  turbulent  rather  than  laminar 
flow.  Nevertheless,  the  injected  air-water  droplet  mixture  is  deflected  by  the  free 
stream  due  to  its  high  velocity  (Miniet=2).  This  is  further  reflected  by  the  volume 
fraction  fields  (Figs.  1 8(b)-2 1  (b)).  As  shown,  the  droplets  carried  by  the  high  velocity 
stream,  spread  as  they  move  towards  the  domain  outlet.  The  gas  temperature  (Fig. 
1 8(c)-2 1  (c))  increases  in  the  region  around  the  injector  due  to  the  large  decrease  in 
velocity  of  the  gas  phase.  This  sudden  decrease  is  supposed  to  be  associated  with  a 
shock  wave.  Evidence  of  this  shock  wave  can  be  noticed  from  the  pressure  contours 
(Figs.  1 8(e)-2 1  (e)),  which  also  shows  the  reflection  of  this  wave  after  hitting  the 
opposite  wall.  In  fact  the  Mach  number  in  the  region  around  the  inlet  decreases  to  very 
low  values.  It  is  suspected  that  smearing  of  this  shock  wave  is  caused  by  numerical 
diffusion,  which  is  a  characteristic  of  the  upwind  scheme.  Unfortunately,  it  has  not 
been  possible  so  far  to  obtain  a  converged  solution  using  any  of  the  high  resolution 
schemes;  this  is  expected  to  be  resolved  with  more  numerical  experimentation. 
However  for  development  purposes,  it  was  deemed  sufficient  to  obtain  results  with  a 
first  order  scheme.  Future  developments  will  address  this  numerical  issue.  The 
predictions  are  physically  sound  and  reveal  a  decrease  in  temperature  due  to 
evaporation  in  regions  where  droplets  exist  (Fig.  1 8(c)-2 1  (c))  and  an  increase  in  the 
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mass  fraction  of  the  water  vapor  in  the  gas  phase  (Fig.  1 8(d)-2 1  (d))  with  increasing 
distance  from  the  injector. 

The  increase  in  the  droplet  temperature  and  decrease  in  the  droplet  diameter  along  the 
centerline  of  the  channel,  as  predicted  with  the  k-e  model,  are  depicted  in  Fig.  22.  It 
should  be  mentioned  here  that  droplets  are  injected  10  cm  from  the  channel  inlet.  Due 
to  the  high  speed  of  the  flow,  the  droplet  temperature  does  not  reach  its  limiting  wet 
bulb  value.  Moreover,  computations  show  that  only  about  5%  of  the  water  droplets 
entering  the  domain  are  transfonned  into  water  vapor  due  to  the  high  velocities 
involved  and  consequently  the  short  residence  time  of  water  droplets. 


Injection  at  Right  Angle  From  Two  Opposing 
Nozzles 


The  physical  situations  for  this  problem,  depicted  in  Figs.  4(c),  is  similar  to  that  of  the 
previous  one,  however  water  is  injected  from  two  nozzles  on  the  lower  and  upper 
horizontal  walls.  Again  the  length  of  the  domain  is  L  (L=l  m)  and  its  width  is  W 
(W=0.25m).  Moreover,  the  required  temperature  dependent  properties  of  the  carrier 
gas  (air)  and  vapor  species  are  given  by  Eq.(121). 

The  Mach  number  and  temperature  of  the  air  at  inlet  to  the  domain  are  taken  to  be  2 
(Mair!iniet=2)  and  900  K,  respectively.  The  mixture  of  air  and  droplet  enters  the  domain 
with  a  temperature  of  300  °K  with  the  volume  fraction  of  the  water  in  the  injected  air- 
water  mixture  being  5x1 0'3,  the  velocity  of  the  fuel-air  mixture  is  10%  the  gas  velocity 
in  the  main  gas  stream,  with  a  fuel  mass  flow  rate  of  2.398  Kg/s.  An  81x41  grid 
system  is  used  and,  as  in  the  previous  case,  the  fuel-air  mixture  is  injected  vertically 
upward  through  a  2mm  nozzle  (subdivided  into  5  control  volumes)  along  the  lower 
horizontal  wall  and  vertically  downward  through  another  2mm  nozzle  (also 
subdivided  into  5  control  volumes)  along  the  upper  horizontal  wall,  with  the  two  jets 
directly  facing  each  others.  Again,  results  are  generated  for  laminar  and  turbulent  flow 
using  the  k-e,  k-co,  and  SST  turbulence  models  and  are  displayed  in  Figs.  23-27. 

Results  follow  the  same  trend  noticed  earlier.  The  droplet  velocity  Figs.  23(a)-26(a) 
and  volume  fraction  fields  depicted  in  Figs.  23(b)-26(b)  reveal  that  both  air-water 
droplet  jets  are  deflected  by  the  free  stream  and  interact  downstream  towards  the  exit 
from  the  channel,  with  the  two  streams  merging  into  one  in  the  laminar  case  only  (Fig. 
23(a))  due  to  the  reasons  stated  above.  The  gas  temperature  increases  in  the  region 
around  the  injectors  due  to  the  formation  of  two  shock  waves  which  are  smeared  by 
the  numerical  diffusion,  and  then  decreases  in  the  area  behind  the  shocks  due  to  the 
heat  absorbed  for  evaporation  (Fig.  23(c)-26(c)).  Moreover,  viscous  dissipation  results 
in  higher  temperatures  in  the  vicinity  of  the  walls.  The  water  vapor  mass  fraction 
increases  as  the  mixture  moves  downstream  in  the  channel  (Figs.  23(d)-26(d))  due  to 
the  increase  in  the  evaporated  amount  with  distance.  Pressure  contours  depicted  in 
Figs.  23(e)-26(e)  reflect  the  formation  of  the  two  oblique  shock  waves  that  are 
reflected  back  and  forth  along  the  upper  and  lower  walls.  Since  these  shock  waves 
move  in  opposite  direction  they  result  in  the  presented  pressure  map. 

The  increase  in  the  droplet  temperature  and  decrease  in  the  droplet  diameter  along  the 
centerline  of  the  channel,  as  predicted  with  the  k-8  model,  are  depicted  in  Fig.  27. 
Again  it  should  be  mentioned  that  droplets  are  injected  10  cm  from  the  channel  inlet. 
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Closing  Remarks 

An  Eulerian  model  involving  discrete  and  continuous  phases  for  the  simulation  of 
droplet  evaporation  and  mixing  at  all  speed  was  fonnulated  and  implemented.  The 
model  allows  for  continuous  droplet  size  changes  without  recourse  to  a  stochastic 
approach.  The  numerical  procedures  are  based  on  a  pressure-based  multi-fluid  finite 
volume  method  and  form  a  solid  base  for  the  future  inclusion  of  combustion  and  other 
modes  of  interactions  such  as  droplet  coagulation  and  break-up.  In  addition  to  the  k-e 
and  k-co  turbulence  models,  the  more  recent  SST  model  was  also  implemented  for  the 
continuous  gas  phase  with  a  coupled  algebraic  model  for  the  discrete  phase.  All 
models  were  appropriately  modified  to  account  for  gas  compressibility  at  high  speeds. 
The  code  was  tested  by  solving  for  evaporation  and  mixing  in  three  very  tough 
configurations,  in  the  various  Mach  number  regimes,  and  over  a  wide  range  of 
operating  conditions  for  both  laminar  and  turbulent  flows.  Even  though  water  droplets 
were  selected  for  the  liquid  phase,  and  only  a  two-component  gas  phase  was  used,  this 
is  not  a  code  or  a  model  limitation.  The  model  can  be  used  with  more  complex  liquid 
fuels,  and  the  code  can  handle  multi-component  multi-fluid  flows,  and  as  such  set  the 
stage  for  the  simulation  of  scramjet  combustion. 
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Figure  Captions 

Fig.  1  Inter-penetrating  media. 

Fig.  2  Terminology  of  the  phase  interface:  A  small  area  of  the  evaporative  surface. 

Fig.  3  Uniform  Droplet  Temperature  Model. 

Fig.  4  Physical  domains. 

Fig.  5  Comparison  of  the  single-phase  and  two-phase  flow  (a)  U  gas  velocity,  (b)  gas  temperature,  (c)  gas  density,  and  (d)  gas  turbulent  kinetic  energy 
profiles  at  three  axial  stations  generated  using  109x160  control  volumes. 

Fig.  6  Comparison  of  the  (a)  gas  velocity,  (b)  gas  temperature,  (c)  vapor  mass  fraction,  and  (d)  droplet  temperature  profiles  at  three  axial  stations  generated 
using  109x160  and  89x120  control  volumes. 

Fig.  7  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass  fraction  of  water  vapor,  and  (e)  pressure  distribution  over  the  domain 

(laminar  flow,  M;n=0.2,  Tgas,m=900  K,  Td,in=300  K,  and  ocd=0.01). 

Fig.  8  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass  fraction  of  water  vapor,  and  (e)  pressure  distribution  over  the  domain 

(turbulent  flow,  k-e  model,  Min=0.2,  Tgasjn=900  K,  Td,in=300  K,  and  ad=0.01). 

Fig.  9  Variation  of  droplet  temperature  and  diameter  along  the  center  line  of  the  domain  (turbulent  flow,  k-e  model,  Min=0.2,  Tgas  ]n=900  K,  Td,in=300  K, 
and  ocd=0.01). 

Fig.  10  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass 

(laminar  flow,  Min=2,  TgaS  ]n=900  K,  Td,in=300  K,  and  oCd=0.01). 

Fig.  11  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass 
(turbulent  flow,  k-e  model,  Min=2,  Tgas  ]n=900  K,  Td,m=300  K,  and  ad=0.01). 

Fig.  12  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass 
(turbulent  flow,  k-co  model,  Mm=2,  Tgasjn=900  K,  Td,in=300  K,  and  ocd=0.01). 

Fig.  13  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass 
(turbulent  flow,  SST  model,  Mm=2,  Tgasm=900  K,  Td,in=300  K,  and  ocd=0.01). 

Fig.  14  Variation  of  droplet  temperature  and  diameter  along  the  center  line  of  the  domain  (turbulent  flow,  k-e  model,  Mm=2,  Tgasjn=900  K,  Td,in=300  K,  and 
ocd=0.01). 

Fig.  15  The  percentage,  by  mass,  of  the  injected  liquid  droplets  that  evaporates  at  different  inlet  gas  temperatures(turbulent  flow,  SST  model,  Mm=2, 
Td,in=300  K,  and  ad=0.01). 
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Fig.  16  The  percentage,  by  mass,  of  the  injected  liquid  droplets  that  evaporates  at  different  inlet  droplet  temperature  (turbulent  flow,  SST  model,  Mm=2, 
Tgas,in=900  K,  and  ad=0.01). 

Fig.  17  The  percentage,  by  mass,  of  the  injected  liquid  droplets  that  evaporates  for  different  domain  length  (turbulent  flow,  k-e  model,  Min=2,  TgaS  ]n=900  K, 
Td,in=370  K,  and  ocd=0.01). 

Fig.  18  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass  fraction  of  water  vapor,  and  (e)  pressure  distribution  over  the  domain 
(laminar  Min=2,  Tgasjn=900  K,  Td,in=300  K,  and  ad=0.01). 

Fig.  19  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass  fraction  of  water  vapor,  and  (e)  pressure  distribution  over  the  domain 
(turbulent  flow,  k-e  model,  Mm=2,  Tgas,in=900  K,  Td,n=300  K,  and  ad=0.01). 

Fig.  20  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass  fraction  of  water  vapor,  and  (e)  pressure  distribution  over  the  domain 
(turbulent  flow,  k-co  model,  Min=2,  Tgasjn=900  K,  Td,in=300  K,  and  ad=0.01). 

Fig.  21  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass  fraction  of  water  vapor,  and  (e)  pressure  distribution  over  the  domain 
(turbulent  flow,  SST  model,  Min=2,  Tgas  in=900  K,  Td;in=300  K,  and  ad=0.01). 

Fig.  22  Variation  of  droplet  temperature  and  diameter  along  the  center  line  of  the  domain  (turbulent  flow,  k-e  model,  Min=2,  TgaS  ]n=900  K,  Td  in=300  K,  and 
ocd=0.01). 

Fig.  23  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass  fraction  of  water  vapor,  and  (e)  pressure  distribution  over  the  domain 
(laminar,  Min=2,  TgaS  ]n=900  K,  Td  m=300  K,  and  ad=0.005). 

Fig.  24  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass  fraction  of  water  vapor,  and  (e)  pressure  distribution  over  the  domain 
(turbulent  flow,  k-e  model,  Min=2,  TgaS  ]n=900  K,  Td  m=300  K,  and  ad=0.005). 

Fig.  25  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass  fraction  of  water  vapor,  and  (e)  pressure  distribution  over  the  domain 
(turbulent  flow,  k-co  model,  Mm=2,  Tgasjn=900  K,  Td,jn=300  K,  and  ad=0.005). 

Fig.  26  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass  fraction  of  water  vapor,  and  (e)  pressure  distribution  over  the  domain 
(turbulent  flow,  SST  model,  Mm=2,  Tgasjn=900  K,  Tdjn=300  K,  and  ad=0.005). 

Fig.  27  Variation  of  droplet  temperature  and  diameter  along  the  center  line  of  the  lower  half  of  the  domain  (turbulent  flow,  k-e  model,  M;n=2,  Tgasdn=900  K, 
Td,in=300  K,  and  ad=0.005). 
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Vg  =  volume  of  gas  phase 
Vd  =  volume  of  droplet  phase 
V  =  Vg  +Vd  =  total  volume 

«d  =  Vg/V 


Figure  1 :  interpenetrating  media 
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Figure  2:  Terminology  of  the  phase  interface 
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Figure  3:  Uniform  Droplet  Temperature  Model. 


Supersonic  Turbulent  Mixing  and  Evaporation 


1 


(a) 


(b) 


Gas  Density 


Gas  Turbulent  Kinetic  Energy 


Figure  5  Comparison  of  the  single-phase  and  two-phase  flow  (a)  U  gas  velocity,  (b)  gas  temperature,  (c)  gas  density,  and  (d)  gas  turbulent  kinetic  energy 
profiles  at  three  axial  stations  generated  using  109x160  control  volumes. 
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Figure  6  Comparison  of  the  (a)  gas  velocity,  (b)  gas  temperature,  (c)  vapor  mass  fraction,  and  (d)  droplet  temperature  profiles  at  three  axial  stations 
generated  using  109x160  and  89x120  control  volumes. 
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Fig.  7  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass 
fraction  of  water  vapor,  and  (e)  pressure  distribution  over  the  domain  (laminar  flow, 
Min=0.2,  TgaSjin=900  K,  Tdiin=300  K,  and  ad=0.01). 
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Fig.  8  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass 
fraction  of  water  vapor,  and  (e)  pressure  distribution  over  the  domain  (turbulent  flow, 
k-e  model,  Min=0.2,  TgaSjin=900  K,  Tdiin=300  K,  and  ad=0.01). 
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Fig.  9  Variation  of  droplet  temperature  and  diameter  along  the  center  line  of  the 
domain  (turbulent  flow,  k-e  model,  Min=0.2,  Tgas  ]n=900  K,  Td,in=300  K,  and  ad=0.01). 


Fig.  10  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass 
fraction  of  water  vapor,  and  (e)  pressure  distribution  over  the  domain  (laminar  flow, 
Min=2,  Tgas  in=900  K,  Tdjin=300  K,  and  (Xd=0.01). 


Fig.  1 1  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass 
fraction  of  water  vapor,  and  (e)  pressure  distribution  over  the  domain  (turbulent  flow, 
k-8  model,  Min=2,  TgaS;in=900  K,  Td,in=300  K,  and  ocd=0.01). 


Fig.  12  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass 
fraction  of  water  vapor,  and  (e)  pressure  distribution  over  the  domain  (turbulent  flow, 
k-co  model,  M,n=2,  TgaS:,n=900  K,  Td,in=300  K,  and  ad=0.01). 


Fig.  13  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass 
fraction  of  water  vapor,  and  (e)  pressure  distribution  over  the  domain  (turbulent  flow, 
SST  model,  Min=2,  Tgas,in=900  K,  Tdjin=300  K,  and  ad=0.01). 
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Fig.  14  Variation  of  droplet  temperature  and  diameter  along  the  center  line  of  the 
domain  (turbulent  flow,  k-e  model,  Min=2,  Tgasjn=900  K,  Td,in=300  K,  and  ad=0.01). 
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Fig.  15  The  percentage,  by  mass,  of  the  injected  liquid  droplets  that  evaporates  at 
different  inlet  gas  temperatures(turbulent  flow,  SST  model,  Min=2,  Td.in=300  K,  and 
ocd=0.01).. 


Fig.  16  The  percentage,  by  mass,  of  the  injected  liquid  droplets  that  evaporates  at 
different  inlet  droplet  temperature  (turbulent  flow,  SST  model,  Mm=2,  Tgas,in=900  K, 
and  ad=0.01).. 
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Fig.  17  The  percentage,  by  mass,  of  the  injected  liquid  droplets  that  evaporates  for 
different  domain  length  (turbulent  flow,  k-e  model,  Min=2,  Tgasjn=900  K,  Tdiin=370  K, 
and  ad=0.01).. 
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Fig.  18  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass 
fraction  of  water  vapor,  and  (e)  pressure  distribution  over  the  domain  (laminar  Mm=2, 
Tgas,in=900  K,  Tdiin=300  K,  and  oCd=0.01). 


Supersonic  Turbulent  Mixing  and  Evaporation 


2 


Fig.  19  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass 
fraction  of  water  vapor,  and  (e)  pressure  distribution  over  the  domain  (turbulent  flow, 
k-8  model,  Min=2,  TgaS,in=900  K,  Td,in=300  K,  and  ocd=0.01). 
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Fig.  20  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass 
fraction  of  water  vapor,  and  (e)  pressure  distribution  over  the  domain  (turbulent  flow, 
k-co  model,  Min=2,  Tgas  ]n=900  K,  Tdiin=300  K,  and  ocd=0.01). 
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Fig.  21  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass 
fraction  of  water  vapor,  and  (e)  pressure  distribution  over  the  domain  (turbulent  flow, 
SST  model,  Mjn=2,  Tgas,in=900  K,  Td,in=300  K,  and  ocd=0.01). 
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Fig.  22  Variation  of  droplet  temperature  and  diameter  along  the  center  line  of  the 
domain  (turbulent  flow,  k-e  model,  Min=2,  Tgasjn=900  K,  Td,in=300  K,  and  ad=0.01). 
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Fig.  23  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass 
fraction  of  water  vapor,  and  (e)  pressure  distribution  over  the  domain  (laminar,  Min=2, 
Tgas,in=900  K,  Td,m=300  K,  and  ad=0.005). 


Supersonic  Turbulent  Mixing  and  Evaporation 


2 


(e) 


Fig.  24  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass 
fraction  of  water  vapor,  and  (e)  pressure  distribution  over  the  domain  (turbulent  flow, 
k-e  model,  Min=2,  TgaS,m=900  K,  Td,in=300  K,  and  ad=0.005). 
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Fig.  25  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass 
fraction  of  water  vapor,  and  (e)  pressure  distribution  over  the  domain  (turbulent  flow, 
k-co  model,  M,n=2,  TgaS:,n=900  K,  Td,in=300  K,  and  ad=0.005). 
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Fig.  26  (a)  droplet  velocity,  (b)  droplet  volume  fraction,  (c)  gas  temperature,  (d)  mass 
fraction  of  water  vapor,  and  (e)  pressure  distribution  over  the  domain  (turbulent  flow, 
SST  model,  Mjn=2,  Tgas,in=900  K,  Td,in=300  K,  and  ad=0.005). 
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Fig.  27  Variation  of  droplet  temperature  and  diameter  along  the  center  line  of  the 
lower  half  of  the  domain  (turbulent  flow,  k-e  model,  Min=2,  Tgas  ]n=900  K,  Tddn=300 

K,  and  ocd=0.005). 


